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ABSTRACT 


Providing  a  simplified  representation  of  terrain  characteristics  has  applications  to 
optimal-path-planning  programs  using  spatial  reasoning.  Utilizing  computer  vision 
techniques,  our  program  creates  polygonal  homogeneous  vegetation  regions  based  on 
map  vegetation  data  from  a  digitized  Defense  Mapping  Agency  database.  Boundary 
points  for  regions  are  identified  from  the  vegetation  codes  in  the  database,  and  then  the 
boundary  contours  of  the  regions  are  traced  using  a  modified  look -left  boundary  tracing 
algorithm.  Each  region  is  then  represented  by  a  polyline  comprised  of  line  segments 
that  meet  a  minimum  threshold  for  fit  using  the  linear  least-squares  criterion.  The 
segments  are  determined  by  first  recursively  splitting  the  region  boundary  until  all 
segments  meet  the  fit  threshold,  and  then  merging  adjacent  segments  that  meet  the 
threshold. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Recent  work  at  the  Naval  Postgraduate  School  has  studied  methods  to  determine 
the  optimal  path  between  two  points  in  terrain.  These  methods  assume  that  the  terrain 
can  be  represented  by  homogeneous  polygonal  regions.  Attributes  of  a  region  can 
include  the  degree  of  vegetation,  type  of  vegetation,  soil  composition,  elevation,  slope, 
orientation  of  the  slope,  and  man-made  physical  features  within  the  area.  Currently  the 
optimal-path-planning  research  is  using  artificial  terrain  data.  A  program  that  will 
take  Defense  Mapping  Agency  digitized  terrain  data  and  create  homogeneous  regions 
would  enable  researchers  to  show  the  real-world  feasibility  of  their  path-planning 
approaches.  In  this  thesis,  we  have  developed  tools  to  create  these  regions  from 
vegetation  data  recorded  at  evenly-spaced  sample  points. 

To  create  two-dimensional  regions  from  vegetation  data  at  evenly-spaced  points 
we  used  two  phases.  The  first  phase  identified  the  boundary  points  of  a  homogeneous 
region,  points  halfway  between  adjacent  sampled  points  of  differing  vegetation.  The 
second  phase  used  a  split-and-merge  method  with  linear-least-squares  constraints  to 
reduce  the  complexity  of  the  region  boundaries  (their  number  of  vertices).  The 
simplified  regions  enable  us  to  represent  the  terrain  characteristics  in  a  clear  and 
concise  manner. 
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B.  ORGANIZATION 


Chapter  2  introduces  previous  work  in  the  areas  of  region  representation  by 
polylines,  split-and-merge  algorithms,  and  the  linear-least-squares  constraints  mefod  of 
approximating  the  fit  of  points  to  a  line.  Chapter  3  describes  the  Defense  Mapping 
Agency  Databases,  and  databases  used  for  optimal  path  planning.  In  Chapter  4  we 
discuss  in  detail  our  program.  Chapter  5  shows  our  experimental  results  in  both 
qualitative  and  quantitative  terms.  Finally,  Chapter  6  summarizes  our  contributions  and 
discusses  some  of  the  possible  areas  for  further  research  based  on  this  work. 
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n.  REVIEW  OF  REGION  REPRESENTATION  BY  POLYGONS 


A.  POLYLINES 

Polylines  can  be  used  to  approximate  the  boundary  of  a  region.  A  polyline 
representation  consists  of  a  list  of  points.  See  Figure  1.  A  region  can  be  represented 


to  any  degree  of  accuracy  by  the  polyline  depending  upon  where  and  how  many 
points  are  used  [Ref.  T.p.  232]. 

A  common  technique  used  to  determine  a  polyline  representing  a  region  is 
splitting  and  merging.  It  can  be  proved  that  the  number  of  line  segments  in  a  polyline 
constructed  using  a  split  and  merge  algorithm  will  never  be  greater  than  two  times  the 
minimum  number  of  segments  for  a  given  line  fit  criteria  [Ref.  2:  p.283].  There  are 
numerous  methods  to  determine  whether  a  given  line  segment  should  be  split,  merged 
or  left  untouched.  Usually  all  methods  split  segments  of  a  polyline  until  all 
subsegments  meet  criteria  for  linearity.  Then  usually  adjacent  segments  of  the  polyline 
are  merged  together  if  they  meet  other  criteria.  The  methodology  used  to  determine 
linearity  is  problem-domain-dependent. 

The  linear  least-squares  method  to  determine  the  collinearity  of  points  has  been 
used  in  a  merge -only  scheme  [Ref.  3].  Two  adjacent  points  were  arbitrarily  selected 
and  points  added  to  their  segment  until  the  linear  least-square  fit  of  the  segment  failed 
the  fit  criteria.  At  that  point  a  new  segment  was  started  and  the  process  repeated  until 
the  region  was  completed.  The  start  point  for  each  segment  determines  the  breakpoints 
of  the  region.  The  primary  disadvantage  of  this  method  is  that  the  merging  without 
backtracking  or  splitting  does  not  give  the  most  aesthetically  pleasing  breakpoints. 

B.  APPROXIMATION  OF  FIT  BY  LINEAR  LEAST-SQUARES  METHOD 

The  fit  of  a  set  of  points  with  respect  to  a  given  line  can  be  judged  using  several 
techniques.  The  linear  least-squares  method  sums  the  squares  of  the  distances  from 
each  data  point  to  the  line.  In  our  work  the  line  is  determined  by  the  two  endpoints 
of  the  boundary-points  subset  tested.  We  use  the  general  equation  of  a  line: 
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Ax  +  By  +  C  =  0 


(1) 


Given  two  endpoints  (XI, Yl)  and  (X2,Y2)  of  the  line  segment  we  can  derive  the  value 
of  the  constants  A,  B,  and  C  by: 

Y  -  Yl  =  Y2  -  Yl  (2) 

X  -  XI  X2  -  XI 

Cross  multiplying  and  solving  for  zero  yields: 


A  =  Yl  -  Y2 

(3) 

B  =  X2  -XI 

(4) 

C  =  X1Y2  -  X2Y 1 

(5) 

The  slope  and  the  intercept  are: 

M  =  -A  /  B 

(6) 

B  =  -C  /  B 

(7) 

The  linear  least-squares  fit  for  n  points  with  respect  to  the  line  calculated  above  is: 


FIT  = 


LU(Ax,^By,  +  Cf 
\  n  (A7  +  B7) 


(8) 
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The  fit  is  a  calculation  of  how  far  on  the  average  the  n  points  are  from  the  line.  A 
fit  of  zero  means  that  all  of  the  points  lie  on  the  line.  By  squaring  the  distance  from 
the  line  we  give  equal  weight  to  points  that  are  on  each  side  of  the  line.  The  square 
root  enables  us  to  renormalize  distance  from  the  line. 
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m.  SPATIAL  TERRAIN  DATABASES 


A.  DEFENSE  MAPPING  AGENCY  DATABASES 

The  Defense  Mapping  Agency  maintains  digitized  terrain  databases  for  most  of 
the  world.  These  databases  contain  information  about  elevation  data,  vegetation,  bodies 
of  water,  and  man-made  objects.  Typically  these  databases  record  data  at  evenly 
spaced  sample  points  in  latitude  and  longitude. 

The  Digital  Terrain  and  Elevation  Data  (DTED)  database  used  in  our  program 
contains  two  types  of  such  infonnation.  The  terrain  elevation  and  the  height  of  the 
vegetation  coverage  are  encoded  in  two  bytes  for  each  sample  point.  The  three  most 
significant  bits  are  the  vegetation  code  and  the  remaining  13  bits  are  the  elevation  in 
feet.  The  three  bits  of  vegetation  code  are  explained  in  Figure  2.  Vegetation  codes 
6  and  7  never  occured  in  the  terrain  samples  we  selected  so  we  did  not  provide 
additional  handling  for  these  codes.  The  samples  are  conducted  every  12.5  meters. 
Each  square  kilometer  or  grid  square  requires  6400  samples  (80  *  80)  or  12,600  bytes. 

B.  DATABASE  USAGE  FOR  OPTIMAL-PATH  PLANNING 

There  are  several  approaches  to  determine  the  optimum  path  between  two  points. 
At  the  Naval  Postgraduate  School  considerable  research  is  being  conducted  using  spatial 
reasoning  Spatial  reasoning  methods  do  not  use  traditional  grid  terrain  modeling 
where  the  terrain  is  represented  by  evenly  distributed  sample  points.  Instead  spatial 
reasoning  uses  descriptive  terrain  modeling  where  the  terrain  is  partitioned  into 
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VEGETATION  CODE 

VEGETATION  HEIGHT 

0 

LESS  THAN  1  METER 

1 

1  -  4  METERS 

2 

4  -  8  METERS 

3 

8  - 12  METERS 

4 

12  -  20  METERS 

5 

GREATER  THAN  20  METERS 

6 

NO  DATA  AVAILABLE 

7 

NOT  USED 

Figure  2  Vegetation  Codes 


homogeneous  regions.  The  minimal-energy  optimal-path-planning  research  conducted 
by  Ron  Ross  uses  partitioned  regions  based  on  the  slope  of  the  terrain,  soil 
composition,  and  other  factors  [Ref  4], 

The  first  step  to  implement  this  work  was  done  by  Seung  Hee  Yee  who  wrote 
a  program  for  planar-patch  terrain  modeling  based  on  the  elevation  data  [Ref  5].  One 
of  the  three  methods  that  he  tested  was  joint  top-down  and  bottom-up  terrain 
modeling.  His  top-down  phase  used  a  quadtree  subdivision  method  to  divide  the 
terrain  area  into  regions  represented  by  a  plane  and  the  fit  of  that  plane  to  the  data 
points.  After  all  subregions  meet  the  appropriate  fit  threshold,  he  used  a  bottom-up 
approach  to  merge  similar  adjacent  regions.  His  bottom-up  phase  used  two  merging 
criteria.  First,  the  adjacent  planes  must  have  a  minimum  difference  of  the  squares  of 
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the  differences  of  the  respective  plane  coefficients.  Second,  the  data  points  from  both 
regions  must  be  representable  by  a  plane  which  meets  the  same  fit  threshold  as  in  the 
top-down  phase.  If  both  criteria  are  met,  then  the  regions  are  merged  and  the  new 
plane  is  stored  along  with  its  fit.  During  the  year  since  Seung  Hee  Yee  developed 
his  program,  Professor  Rowe  has  improved  it  by  including  better  techniques  to  insure 
continuity  between  the  planes  of  adjacent  regions.  The  planar  patches  or  regions  from 
these  programs  could  be  further  partitioned  into  regions  of  homogeneous  vegetation,  but 
this  has  not  yet  been  done. 

The  simplest  minimal-energy  optimal -path -planning  programs  run  on  the  order  of 
N  squared  with  respect  to  the  number  of  vertices  per  region  and  the  number  of  regions. 
For  this  reason  it  is  important  to  create  regions  which  preserve  the  topology  of  the 
original  sampled  data  and  yet  minimize  the  number  of  vertices. 
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IV.  IMPLEMENTATION 


The  program  developed  is  a  method  to  create  polygonal  homogeneous-vegetation 
regions  from  the  Defense  Mapping  Agency  database  described  in  Chapter  3.  To  form 
these  regions  the  program’s  first  phase  identifies  which  sample  points  in  the  database 
have  different  vegetation  codes  than  their  neighbors.  These  sample  points  are  identified 
as  boundary  cells.  A  boundary-identification  algorithm  is  then  used  on  these  boundary 
cells  to  determine  sequences  of  boundary  points  between  regions  of  different  vegetation. 
A  second  phase  then  takes  the  sequence  of  boundary  points  for  each  region  and 
performs  a  split-and-merge  algorithm  on  the  list  of  those  points  to  redefine  the  regions 
in  a  more  efficient  form  (one  with  fewer  vertices). 

A.  DESCRIPTION  OF  THE  BOUNDARY-IDENTIFICATION  PHASE 

We  chose  C  as  the  language  for  the  boundary-identification  phase.  The  program 
runs  using  a  BSD  4.3  compiler  on  a  VAX  11/785  or  on  an  Integrated  Solutions  (ISI) 
workstation.  Appendix  B  contains  the  source  code  for  the  boundary  indentification 
phase.  This  phase  requires  input  from  the  Defense  Mapping  Agency  DTED  database. 
This  database  consists  of  elevation  and  vegetation  data  points  every  12.5  meters.  Our 
program  builds  regions  containing  sample  points  with  the  same  vegetation  code  for  a 
one  kilometer  by  one  kilometer  grid  square. 

Several  assumptions  were  required  in  the  boundary-identification  phase.  First, 
each  region  formed  is  assumed  to  contain  vegetation  of  only  one  kind.  This  is 
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fundamental  to  the  concept  of  homogeneous  regions.  Second,  regions  will  not  consist 
of  one  sample  point,  These  singular  points  which  have  no  adjacent  neighbors  above, 
below,  or  beside  them  are  changed  to  have  die  same  code  as  the  adjacent  data  point 
with  the  largest  vegetation  code  (i.e.  heaviest  vegetation).  Third,  holes  in  regions  will 
be  detected  only  in  the  sense  that  another  region  will  be  formed  inside  of  the  outer 
region.  No  special  classification  is  added  to  a  region  if  it  has  a  hole. 

The  boundary-identification  phase  is  divided  into  two  passes  through  the  input 
data.  The  first  pass  identifies  the  boundary  cells.  A  sample  point  is  considered  to  be 
a  boundary  cell  if  any  of  its  four  neighbors  above,  below,  or  beside  it  have  a  different 
vegetation  code.  The  sole  exception  to  this  rule  exists  when  the  point  itself  or  its 
differing  neighbor  are  singular  points,  which  are  not  considered  to  be  boundary  cells 
in  accordance  with  the  second  assumption.  An  example  of  vegetation  codes  and 
identified  boundary  cells  is  shown  in  Figure  3. 

The  second  pass  searches  through  the  identified  boundary  cells  and  traces  the 
contour  for  each  region:  as  boundary  cells  are  traced,  they  are  marked  and  identified 
with  a  specific  region  number.  Region  numbers  are  generated  as  necessary  beginning 
with  1.  The  search  for  a  boundary  cell  to  start  the  contour-tracing  for  each  region  is 
started  in  the  upper  left-hand  comer  of  the  input.  The  grid  is  searched  in  a  left-to- 
right  raster  scan  stopping  at  the  first  boundary  cell  that  has  not  been  assigned  to  a 
specific  region.  After  this  region  has  been  traced  by  the  contour-tracing  algorithm  from 
this  start  point,  the  raster  scan  continues. 

The  contour-tracing  algorithm  follows  the  traditional  eight-connected  look-left 
algorithm  for  traversing  a  maze  [Ref.  6:p.  278],  At  each  boundary  cell  during  the  trace 
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Figure  3  Boundary  Cells 


the  direction  of  entry  determines  which  way  we  look  to  find  the  next  boundary  cell; 


the  direction  of  entry  at  the  start  point  is  defined  to  be  the  raster  scan  direction. 


Normally,  we  look  to  the  boundary  cell  to  the  left  of  the  direction  of  entry.  If  this 


cell  is  not  a  boundary  cell  of  the  same  vegetation  code  we  look  one  cell  left  and 
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forward.  This  process  is  continued  clockwise  until  we  find  a  boundary  cell.  The 
contour-tracing  algorithm  is  continued  until  we  reach  the  start  point  or  a  border  of  the 
picture  (the  first  pass  ensures  that  either  event  must  eventually  occur). 

Several  additions  to  the  contour-tracing  algorithm  were  included  to  enable  us  to 
handle  conditions  along  the  exterior  borders  of  the  overall  1km.  by  1km.  area.  When 
the  contour  tracing  reaches  an  exterior  border,  a  check  is  made  to  determine  if  a  border 
has  previously  been  encountered  for  the  current  region.  If  a  border  has  not  previously 
been  encountered,  the  tracing  continues  in  the  opposite  direction  starting  at  the  initial 
point;  otherwise,  tracing  stops.  See  Figure  4.  Second,  if  the  initial  point  of  the  region 
is  next  to  the  exterior  border  the  tracing  will  start  away  from  the  border.  See  Figure 
5. 

A  significant  improvement  in  the  representation  of  the  data  during  this  pass  is 
gained  by  treating  the  actual  boundary  between  two  regions  as  the  set  of  all  points 
half-way  between  the  centers  of  each  pair  of  adjoining  boundary  cells  with  a  different 
vegetation  code.  Tracing  thus  proceeds  between  these  points.  This  is  similar  to  the 
crack  edges  approach  used  to  represent  the  boundary  between  regions  in  [Ref.  l:p.  78]. 
Additionally,  this  teclinique  smooths  the  staircase  effect  created  by  the  crack  edges 
approach.  See  Figure  6. 

The  boundary -tracing  phase  trades  memory  space  for  clarity  of  code.  The  input 
is  stored  in  an  80  x  80  array  which  represents  the  vegetation  codes  at  each  of  the 
6400  sampled  points  in  row-major  order.  The  first  pass  places  boundary  cells  in  a 
6400  x  5  array  in  row-major  order,  classifying  each  by  its  x-coordinate,  y-coordinate, 
vegetation  code,  whether  it  is  a  boundary  point,  and  a  flag  for  marking  when  the 


Figure  4  Tracing  after  Hitting  the  Border 


second  pass  places  the  boundary  point  in  a  specific  region.  The  second  pass  stores 
into  a  separate  array  for  each  region  the  sequenced  boundary  points  x -coordinate,  y* 
coordinate,  and  vegetation  code.  The  array  of  regions  is  output  to  a  formatted  file 
which  will  be  used  as  input  to  the  split-and-merge  phase  described  in  section  B. 


FIGURE  5  Conditions  for  Starting  Away  from  the  Border 


The  boundary -identification  phase  is  limited  to  grids  containing  no  more  than 
6400  input  data  points,  and  can  handle  up  to  75  output  regions  of  no  more  than  400 
traced  points.  These  constraints  were  selected  to  insure  the  program  can  handle  grid 
squares  with  either  many  small  regions,  or  grid  squares  with  a  few  large  regions. 
Although  the  normal  cases  for  the  boundary  tracing  algorithm  have  been  tested,  the 


Figure  6  Boundary  Point  Representation  of  Boundary  Between  Regions 


handling  of  all  possible  cases  involving  small  regions,  such  as  those  containing  2 
boundary  cells  near  other  regions,  has  not  been  verified  to  work,  as  in  some  unusual 
situations  the  tracing  could  skip  points. 
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B.  DESCRIPTION  OF  THE  SPLIT-AND-MERGE  PHASE 


The  split-and-merge  phase  of  the  program  is  written  in  the  M-Prolog 
programming  language  for  the  Integrated  Solutions  (ISI)  workstations.  The  list¬ 
processing  capabilities  of  Prolog  enables  us  to  treat  each  region  as  a  list  of  points  and 
conduct  operations  on  the  members  of  the  list. 

The  input  to  the  split-and-merge  phase  is  a  Prolog  fact  for  each  region  detected 
in  the  boundary-identification  phase.  Each  fact  contains  an  identification  number  of 
the  region,  die  list  of  points  for  the  region,  the  vegetation  code  of  the  region,  and  the 
total  number  of  points  in  the  region.  See  Figure  7.  The  output  of  the  split-and-merge 
phase  is  a  set  of  facts  which  contain  vegetation  codes  and  lists  of  the  coordinates  of 
each  vertex  for  the  polygonal  homogeneous  vegetation  regions. 

The  major  data  structures  used  in  the  split-and-merge  phase  are  lists,  expressed 
as  facts  asserted  throughout  the  phase.  A  global  variable  is  created  for  the  processing 
of  each  region.  This  contains  the  list  of  boundary  points  for  that  region;  indexing 
based  on  the  placement  of  the  boundaty  points  in  the  list  is  used  to  access  it  in  the 
linear  least-squares  module.  This  saves  allocated  memory  in  terms  of  the  statement 
table,  global  stack,  and  the  main  stack.  A  segment  fact  is  asserted  each  time  the 
linear  least-squares  Fit  is  calculated  for  a  line  segment  during  the  splitting  process  of 
the  split-and-merge  phase.  See  Figure  9.  These  facts  are  asserted  and  retracted  often 
as  the  program  seeks  the  proper  combination  of  segments. 

The  split-and-merge  phase  contains  five  modules;  the  source  code  is  contained 
in  Appendix  C.  The  split-and-merge  module  controls  the  phase  using  the  built-in 
automatic  backtracking  of  Prolog.  The  linear  least -squares  module  calculates  the  fit 
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module  regions. 

/*$eject*/ 

body. 

reg(0,[[0.0e0,22.5e0], 

[1.0e0,21.0e0], 

[1.0e0,21.5e0], 

[2.0e0,21.5e0], 

• 

• 

[43.5e0,0.0e0]].3,62). 

reg(l,  [[0.0e0,72.5e0], 

• 

• 

• 

/*  The  first  fact  describes  region  number  0. 

*/ 

/*  The  region  contains  62  data  points. 

*1 

/*  The  value  for  the  region  is  vegetation  code  3. 

*/ 

Figure  7  Input  Facts  for  Split-and-Merge  Phase 


of  points  in  the  vicinity  of  a  line.  The  list  module  provides  basic  Prolog  list¬ 
processing  predicates  The  math  module  provides  math  funtions  such  as  the  square 
root  which  are  not  provided  in  M-Prolog.1  Finally,  the  regions  module  provides  the 
data  input  from  the  boundary-identification  phase. 

The  split-and-merge  phase  processes  one  pair  of  adjacent  regions  at  a  time, 
starting  with  pairs  with  non-zero  vegetation  codes.  For  each  adjacent  region,  the  list 
of  boundary  points  from  the  two  adjacent  regions  is  intersected  to  produce  a  list  of 
adjacent  boundary  points.  To  achieve  clean  boundaries  between  the  two  regions,  no 


'  The  source  code  for  the  linear  least-squares  module  and  most  of  the  math  and  list 
modules  were  written  by  Professor  Rowe.  Prolog  predicates  written  by  Professor  Rowe  are 
marked  with  an  asterisk  in  Appendix  C. 
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asserta(segment(R,N  1  ,N2  JFit  1 2» . 


f*  R  represents  the  region  number  of  the  segment.  */ 

f*  N1  and  N2  are  the  vertices  of  the  line  segment.  */ 

/*  N1  and  N2  are  integer  indices  into  the  list  of  points  */ 

/*  for  the  region.  */ 

/*  Fit  12  is  the  linear  least-squares  fit  of  the  ponts  between  N1  and  */ 

f*  and  N2  and  the  line  defined  by  N1  and  N2.  */ 


Figure  8  Segment  Facts 

merging  of  line  segments  is  allowed  to  use  points  from  boundary  point  lists  of  differing 
region  pairs.  See  Figure  8.  The  adjacent  point  list  between  the  two  regions  will  be 
handled  twice  for  the  two  regions  of  the  pair;  however,  the  list  will  be  split  and 
merged  identically  in  each  case. 

The  splitting  process  follows  the  traditional  divide -and-conquer  algorithm  [Ref  2:p. 
282].  The  linear  least-squares  fit  is  calculated  for  a  segment,  and  a  segment  fact  is 
asserted.  If  the  fit  is  less  than  the  splitting  threshold,  the  segment  is  split  into  two 
new  line  segments  which  share  a  common  point,  the  new  "breakpoint"  between  them, 
and  the  old  segment  fact  is  retracted.  This  process  continues  recursively  until  there 
are  no  more  segment  facts  to  be  split. 

The  calculation  of  the  linear  least-squares  fit  is  done  in  the  linear  least-squares 
module  using  equation  (8)  from  Chapter  2.  The  end  points  of  each  line  segment  define 
the  line  that  will  be  tested  for  the  fit.  The  use  of  the  end  points  to  define  the  line 
segment,  rather  than  searching  for  the  "best"  line,  enables  us  to  connect  adjacent  line 
segments  at  exactly  the  breakpoints. 
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•  POINTS  DEFINING  THE  ADJACENT  POINT  LIST 


Figure  9  Adjacent  Regions 

The  merge  process  examines  each  pair  of  line  segments  with  a  common  index 
number  (breakpoint)  in  the  segment  facts  list.  The  segments  are  provisionally  merged 
and  the  fit  is  calculated.  If  the  fit  is  less  than  or  equal  to  the  merge  threshold,  the 
provisional  merge  is  made  permanent.  A  new  segment  fact  is  then  asserted,  and  the 


two  original  segment  facts  are  then  retracted.  This  process  is  continued  until  no  more 
adjacent  line  segments  can  be  merged. 

The  final  process  in  the  linear  least-squares  module  of  the  split-and-merge  phase 
takes  the  list  of  index  numbers  (representing  the  placement  of  the  boundary  points 
within  the  input  list)  and  finds  their  actual  x  and  y-coordinates.  This  shortened  list  is 
then  bound  to  a  newregion  fact  along  with  the  value  of  the  region. 
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V.  RESULTS 


A.  PERFORMANCE  MEASURES  OF  THE  BOUNDARY  -IDENTIFICATION 
PHASE 

In  test  runs,  the  boundary-identification  phase  appears  to  properly  define  the 
boundaries  of  both  convex  and  concave  regions,  including  regions  with  holes.  The 
maximum  error  for  any  point  on  the  region  boundary  is  equal  to  one  half  of  the 
distance  between  the  sample  points,  ignoring  singular  points.  Since  the  input  data 
points  were  evenly  spaced  every  12.5  meters,  the  maximum  error  at  any  point  is  equal 
to  6.25  meters. 

The  phase  is  memory-intensive.  The  program  allocates  196800  bytes  for  the 
three  arrays  that  are  used  to  temporarily  store  information.  This  is  more  than  three 
times  the  number  of  bytes  required  to  process  any  of  the  test  grid  squares.  Since  we 
are  not  trying  to  optimize  the  code,  we  allocated  more  room  than  would  probably  be 
needed.  Extensive  use  of  file  I/O,  and  the  N  squared  nature  of  the  boundary- 
identification  algorithm,  where  N  is  the  length  of  the  square  grid,  slow  the  program 
to  an  average  run  time  of  1  minute. 

B.  PERFORMANCE  MEASURES  OF  THE  SPLIT-AND-MERGE  PHASE 

The  regions  obtained  from  the  split-and-merge  phase  appear  to  represent  the 

input  database  without  significant  error.  Examples  of  input  terrain  and  the  regions 
formed  are  contained  in  Appendix  A.  As  expected,  the  fit  threshold  determines  the 
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number  of  vertices  for  the  output  regions.  Table  1  shows  several  thresholds  and  the 
resulting  number  of  vertices  for  each  region. 


TABLE  1  EFFECT  OF  VARYING  THE  THRESHOLD  ON  THE  NUMBER  OF 

VERTICES  PER  REGION 


NUMBER  OF  POLYGON  VERTICES 

REGION  # 

THRESHOLD 

THRESHOLD 

THRESHOLD 

THRESHOLD 

0.10 

1.0 

2.50 

10.0 

0 

21 

11 

7 

3 

1 

3 

3 

2 

2 

2 

5 

4 

3 

2 

3 

23 

15 

9 

5 

4 

21 

11 

6 

4 

5 

21 

12 

4 

4 

6 

18 

12 

7 

4 

A  threshold  of  0.01  allows  almost  no  merging.  Conversely,  a  threshold  of  10.0 
reduces  the  regions  to  only  a  few  vertices. 

The  split-and-merge  phase  is  expensive  in  both  time  and  space.  Table  2  shows 
the  maximum  observed  quantities  for  the  Main  Stack,  Global  Stack,  Statement  Table, 
Evaluations,  cpu  run  time,  and  the  real  run  time.  The  stacks  and  statement  table  are 
managed  dynamically  in  our  program,  so  the  numbers  vary  through  the  running  of  the 
program.  The  higher  numbers  were  obtained  while  processing  large  lists  in  the  linear 


TABLE  2  MAXIMUM  OBSERVED  SIZE  FOR  M-PROLOG  SYSTEM 
PARAMETERS  IN  THE  SPLIT-AND-MERGE  PHASE 


PARAMETERS 

OUANTITY 

MAIN  STACK 

1057  ITEMS 

GLOBAL  STACK 

4544  ITEMS 

STATEMENT  TABLE 

66,939  STATEMENTS 

CPU  RUN  TIME 

20  MINUTES 

23  SECONDS 

REAL  RUN  TIME 

30  MINUTES 

least-squares  module.  Prior  to  entering  the  linear  least-squares  module  the  Main 
Stack,  Global  Stack,  and  the  Statement  Table  were  860  items,  660  items,  and  60,836 
statements  respectively.  These  numbers  only  reflect  the  points  in  the  program  that  we 
observed  the  system  parameters.  It  is  possible  that  these  parameters  could  exceed  the 
figures  stated  in  Table  2  in  other  portions  of  the  program. 
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VI.  CONCLUSION 


The  primary  objective  of  this  work  was  twofold.  First,  we  needed  to  show  the 
feasibility  of  creating  a  polygonal  region  representation  of  terrain  vegetation  from 
digitized  map  data.  Our  program  shows  that  this  can  be  done.  The  program  provides 
a  suboptimal  solution  that  visually  appears  to  properly  represent  the  terrain.  But  the 
split-and-merge  phase  of  the  program  is  time-consuming  as  the  program  searches  for 
proper  line  segments.  The  many  calculations  involved  implementing  the  fit  calculations 
in  M-Prolog  is  the  main  reason;  other  less  mathematical  methods  for  determining  the 
fit  of  lines  do  exist  and  could  easily  replace  the  linear  least-squares  calculations  without 
changing  the  rest  of  the  split-and-merge  phase. 

The  second  objective  of  this  work  was  to  provide  a  working  tool  for  path¬ 
planning  research  at  the  Naval  Postgraduate  School.  This  program  will  enable 
researchers  to  create  polygonal  regions  of  terrain  vegetation. 

The  program  only  considers  a  1km.  by  1km.  grid  square.  A  more  useful  program 
for  optimal  path  planning  could  consider  much  larger  areas  of  terrain.  This  would 
probably  require  more  efficient  data  storage  in  the  boundary-identification  phase  and 
possibly  optimizing  and  compiling  the  split-and-merge  phase.  An  alternative  and  less 
costly  technique  would  take  the  results  from  several  adjacent  grids  and  attempt  to 
splice  regions  together  that  hit  their  respective  borders. 

Another  topic  for  future  research  would  be  overlaying  the  results  of  this  work 
with  the  three-dimensional  planar-patch  terrain  models  that  Seung  Hee  Yee  and 
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Professor  Rowe  developed.  This  would  create  partitioned  homogeneous  regions  of 
constant  slope  as  well  as  vegetation. 
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APPENDIX  A  TEST  RESULTS 
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Figure  5.4  Raw  Data  for  35  57’  30  "N,  121  17’  W 
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THRESHOLD  EQUALS  0.01 


0 


NUMBERS  ARE  VEGETATION  CODES 


THRESHOLD  EQUALS  2.50 


0 


NUMBERS  ARE  VEGETATION  CODES 


THRESHOLD  EQUALS  10.0 


NUMBERS  ARE  VEGETATION  CODES 


THRESHOLDS  1.0  AND  0.01  SUPERIMPOSED 
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THRESHOLDS  0.01  AND  2.50  SUPERIMPOSED 
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THRESHOLDS  10.0  AND  0.01  SUPERIMPOSED 


APPENDIX  B  SOURCE  CODE  FOR  THE  BOUNDARY-IDENTIFICATION 

PHASE 


#include  <ctype.h> 

#include  <stdio.h> 

#includc  <math.h> 

#include  "header.h" 

unsigned  short  vegarray[80][80J; 

unsigned  short  bndpts[6400][5];  /*  6400  pts  x,  y,  value,  boundary?,  reg?  */ 

float  regs[75][400][3];  /*  75  regions  x  400  pts  with  xl,  yl  in  array  coords  */ 

/*  and  the  region  vdue.  */ 


unsigned  short  pnt; 
unsigned  short  veg; 
unsigned  short  new_val; 
unsigned  short  curr_dir; 
unsigned  short  curr_pt; 
unsigned  short  curr_x; 
unsigned  short  curr_y; 
unsigned  short  neighbor_x; 
unsigned  short  neighbor_y; 
unsigned  short  ij,k; 
boolean  found_next; 
boolean  hit_border; 
boolean  stop_tracing; 
boolean  trace_counter; 
boolean  trace_border; 
boolean  start_counter; 
int  fdin; 

main() 

I 

unsigned  short  x,y; 
int  is_not_singlept(); 
int  single_adjacen?  pt(); 
int  is_boundarypuy; 
int  can_extend(); 
int  adjacent_border(); 
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int  curr_dir_diagonal(); 
int  start_backwards(); 
fdin  =  open(infile,0); 

/*  creates  the  vegatation  array  80  x  80  with  values  */ 
for  (x=0;  x  <  80;  x-H-){ 
for  (y  =  0;  y  <  80;  y++)( 
read(fdin,  &veg,  2); 
vegarray[x][y]  =  veg  -  48; 

) 

) 

pnt  =  0; 

i  =  0;  /*  The  current  point  in  the  second  pass  */ 
j  =  0;  /*  The  first  four  border  regions  are  initiated  w/out  j  */ 
k  =  0;  /*  Next  point  #  for  current  region  being  traced  */ 
for  (x  =  0;  x  <  80;  x++){ 
for  (y  =  0;  y  <  80;  y++){ 

if  (x  ==  0  II  x  =  79  II  y  =  0  II  y  =  79)  { 
bndpts(pntl[0]  =  x; 
bndpts[pnt][l]  =  y; 
bndpts[pnt][2]  =  vegarray[x][y]; 
bndpts[pnt][3]  =  FALSE; 
bndpts[pnt][4]  =  TRUE; 

} 

else  if  (((vegarray[x][y]  !=  vegarray[x][y+l]  && 
is_not_singlept(x,y+l))  II 
(vegarray[x][y]  !=  vegarray[x+l][y]  && 
is_not_singlept(x+ 1  ,y ) )  II 
(vegarray[x][y]  !=  vegarray[x][y-l]  && 
is_not_singlept(x,y- 1 ))  II 
(vegarray[x][y]  !=  vegarray[x-l][y]  && 
is_not_singlept(x-l,y)))  && 
is_not_s  inglept  (x  ,y ) )  ( 
bndpts[pnt][0J  =  x; 
bndpts[pntj[l  1  =  y; 
bndpts[pnt][2J  =  vegarray[x][y]; 
bndpts[pnt][3]  =  TRUE; 
bndpts[pnt][4]  =  FALSE; 

}  /*  if*/ 

else  if  (single_adjacent_pt(x,y))  { 
bndpts[pnt][0]  =  x; 
bndpts[pntj[l]  =  y; 
bndpts[pnt][2]  =  new_val; 
bndpts[pnt][3]  =  TRUE; 
bndpts[pnt][4]  =  FALSE; 

}  /*  if  */ 
else  ( 
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bndpts[pnt][0]  =  x; 
bndpts[pnt][lj  =  y; 
bndpts[pntj[2]  =  vegarray[x]ly]; 
bndpts[pnt][3]  =  FALSE; 
bndpts[pntj[4]  =  FALSE; 

} 

pnt++; 

j  [*  for  y  */ 

}  /*  for  x  */ 

t*  The  second  Pass  through  the  data  base  starts  here.  This  pass  traces  */ 
f*  the  boundaries  of  each  region  by  scanning  through  the  entire  DB  looking  */ 
/*  for  a  boundary  point.  Once  a  boundary  point  is  found  it  starts  tracing  */ 
f*  around  the  region.  The  output  is  to  the  regs  set  of  arrays.  */ 

while  (i  <  6400)  { 
cuir_dir  =  North; 
k  =  0; 

hit_border  =  FALSE; 
trace_border  =  FALSE; 
if  (is_boundarypt())  { 
curr_pt  =  i; 

curr_x  =  bndpts[curr_pt][0]; 
curr_y  =  bndpts[curr_ptJll]; 
stop_tracing  =  FALSE; 
trace_counter  =  FALSE; 
start_counter  =  FALSE; 
if  (adjacent_border()) 
connect_border(); 
if  (stan_backwards())  { 
start_counter  =  TRUE; 
trace_counter  =  TRUE; 
hit_border  =  TRUE; 
curr_dir  =  SEast; 

I 

do  { 

if  (start_counter)  { 

right_front_neighbor( ); 
add_point(); 
found_next  =  FALSE; 
while  (!found_next)  ( 
if  (can_extend())  { 
move_curr_dir(); 
found_next  =  TRUE; 

I 

else  { 

look_left_one(): 
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if  (can_extend())  { 
move_curr_dir(); 
found_next  =  TRUE; 

} 

else  ( 

front_neighbor(); 
add _point(); 
look_left_one(); 

}  l*  inner  else  */ 

}  /*  outer  else  */ 

}  /*  while  not  found  next  */ 
start_counter  =  FALSE; 

}  /*  if  start  counter  */ 

else  if  (trace_counter  &&  curr_dir_diagonal())  { 
right_rear_neighbor( ); 
add_point(); 
look_right_two(); 
found_next  =  FALSE; 
while  (!found_next){ 
if  (can_extend())  j 
move_curr_dir( ) ; 
found_next  =  TRUE; 

} 

else  { 

look_left_one(); 
if  (can_extend())  { 
move_curr_dir(); 
found_next  =  TRUE; 

) 

else  { 

front_neighbor(); 

add_pointO; 
look_left_one( ): 

)  /*  inner  else  */ 

}  /*  outer  else  */ 

)  /*  while  not  found_next  */ 

if  (curr_pt  ==  i) 
continue; 

)  /*  if  curr  dir  is  diagonal  and  tracing  counter  clockwise  */ 
else  if  (trace_counter  &&  !  curr_dir_diagonal()  && 
!trace_border)  { 
right_neighbor(); 
add_point(); 
look_right_one(); 
found_next  =  FALSE; 
while  (!found_next)  | 
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if  (can_extend())  { 
move_curr_dir( ); 
found_next  =  TRUE; 

) 

else  { 

look_left_one(); 
if  (can_extend())  { 
move_curr_dir(); 
found_next  =  TRUE; 

} 

else  { 

front_neighbor(); 

add_poim(); 

look_left_one(); 

}  /*  inner  else  */ 

)  /*  outer  else  */ 

)  /*  while  not  found_next  */ 
if  (curr_pt  ==  i) 
continue; 

)  /*  else  curr_dir  is  not  diagonal  */ 
else  if  (curr_dir_diagonal())  { 
left_rear_neighbor( ); 
add_point(); 
look_left_two(); 
found_next  =  FALSE; 
while  (!found_next){ 
if  (can_extend())  { 
move_curr_dir(); 
found_next  =  TRUE; 

) 

else  { 

look_right_one(); 
if  (can_extend())  ( 
move_curr_dir( ); 
found_next  =  TRUE; 

) 

else  ( 

front_neighbor(); 

add_point(); 

look_right_one(); 

}/*  inner  else  */ 

)  /*  outer  else  */ 

}  /*  while  not  found_next  */ 

if  (curr_pt  ==  i) 
continue; 

)  /*  if  curr  dir  is  diagonal  */ 


( *  ***********************  eise  curr_dir  is  not  diagonal  ***********  */ 
else  if  (!trace_border)  ( 
left_neighbor(); 
add_point(); 
look_left_one(); 
found_next  =  FALSE; 
while  (!found_next)  { 
if  (can_extend())  { 
move_curr_dir( ) ; 
found_next  =  TRUE; 

} 

else  { 

look_right_oneO ; 
if  (can_extend())  { 
move_curr_dir( ); 
found_next  =  TRUE; 

} 

else  { 

front_neighbor(); 

add_point(); 

look_right_one(); 

}  /*  inner  else  */ 

)  /*  outer  else  */ 

}  /*  while  not  found_next  */ 
if  (curr_pt  ==  i) 
continue; 

)  /*  else  curr_dir  is  not  diagonal  */ 
if  (adjacent_border())  { 

bndpts[curr_pt][4]  =  TRUE; 
connect_border( ); 

)  /*  if  current  point  is  adjacent  border  */ 

}  while  ((!(curr_pt  ==  i  &&  !(trace_counter)))  && 
!(stop_tracing)); 

/ *  If  the  region  has  the  same  stan  and  end  point  make  the  last  point  the  */ 
/*  same  as  the  first  point.  */ 

if  ((curr_pt  ==  i)  &&  (!hit_border)) 
duplicate_first_pt( ); 

/ '*  Set  a  flag  at  the  end  of  each  region  so  the  output  routine  knows  to  stop  */ 
regs(j][k][0]  =  9999; 
regs[j]lk][l]  =  9999; 
j++;  f*  increment  region  number  */ 

}  /*  if  is  boundary  point  */ 

i++;  /*  increment  point  number  that  we  are  checking  for  bndy  */ 

}  /*  while  i  <  6400  */ 
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File  BNDUTIL.C 

/*  This  file  contains  the  necessaiy  functions  to  execute  the  bndpts.c  program.*/ 

♦include  <ctype.h> 

♦include  <stdio.h> 

♦include  "header.h" 

extern  unsigned  short  veg; 

extern  unsigned  short  new_val; 

extern  unsigned  short  curr_dir; 

extern  unsigned  short  curr_pt; 

extern  unsigned  short  curr_x; 

extern  unsigned  short  curr_y; 

extern  unsigned  short  ij,k; 

extern  unsigned  short  neighbor_x; 

extern  unsigned  short  neighbor_y; 

extern  boolean  trace_counter; 

extern  boolean  trace_border; 

extern  boolean  stop_tracing; 

extern  boolean  hit_border; 

extern  unsigned  short  vegarray[80][80]; 

extern  unsigned  short  bndpts[6400][5]; 

extern  float  regs[75][400][3]; 

extern  int  fdout; 

extern  int  is_boundarypt(); 

extern  int  is_singlept(); 

extern  int  single_adjacent_pt(); 

j*p  ******************************************************************/ 

int  is_not_singlept(x,y) 
unsigned  short  x.y: 

{ 

return  (vegarray[x][y]  ==  vegarray[x+l]ly]  II 
vegarray[x][y]  ==  vegarray[x-l][y]  II 
vegarray[x][y]  ==  vegarray[x][y+l]  II 
vegarray[x]ly]  ==  vegarray[x][y-l]); 

}  /*  function  is_not_singlept  */ 

!+  ********************************************************************* 

V 

int  single_adjacent_pt(x,y) 
unsigned  short  x,y; 

I 

/*  First  find  the  maximum  4  way  neighbor  of  the  point.  */ 

newval  =  0; 
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if  (vegarray[x  +  1]  [y]  >  vegarray[x][y  +  1]) 
new_val  =  vegarray[x+l][y]; 
if  (vegarray[x][y-l]  >  vegarray[x+l][y]) 
new_val  =  vegarray[x][y-l]; 
if  (vegarray[x-l][y]  >  vegarray[x][y-l]) 
new_val  =  vegarray[x-l][yj; 
if  (vegarray[x][y+l]  >  vcgarray[x-l][y]) 
new_val  =  vegarray[x][y+lj; 

/*  Second,  return  true  if  the  point  has  4  way  neighbors  with  at  least  two  */ 

/*  different  values.  Returns  largest  neighbors  value  as  new_val.  */ 

retum(!  is_not_singlept(x,y)  && 

(vegarray[x+l][y]  !=  vegarray[x][y+l]  II 
vegarray[x+l][y]  !=  vegarray[x][y-l]  II 
vegarray[x+l][y]  !=  vegarray[x-l]ly])); 

)  /*  function  single_adjacent_pt  */ 

I *  *********** **************** ****************************************** 

V 

/*  A  point  is  a  boundary  point  if  it  is  a  boundary  and  it  is  not  already  */ 

/*  part  of  a  region  and  it  does  not  have  a  value  of  zero.  */ 

int  is_boundarypt() 

I 

return  (bndptsfi][3]  ==  TRUE  &&  bndpts[i][4]  ==  FALSE  && 
bndpts[i][2]  !=  0); 

}  /*  function  is_boundarypt  */ 

^*  ***********4,^********«t********************************************** 

*/ 

/*  A  region  is  extended  if  the  point  in  the  curr  direction  is  a  boundary  */ 

/*  and  the  point  in  the  curr  direction  is  the  same  value  as  the  current  */ 

/*  point.  */ 

v 

int  can_extend() 

I 

switch  (curr_dir) 

I 

case  North,  return  (bndptsfcurr_pt  +  1][3]  =  TRUE  && 

bndpts[curr_pt  +  1][2]  =  bndpts[curr_pt][2]); 
case  NEast:  return  (bndpts[curr_pt  +  81][3]  =  TRUE  && 

bndpts[curr_pt  +  81J[2J  =  bndpts(curr_pt][2]); 
case  East.  return  (bndpts[curr_pt  +  80][3]  =  TTIUE  && 

bndpts[curr_pt  +  80J[2J  ==  bndpts[curr_pt][2]); 
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if  (bndpts[curr_pt][0]  =  1)  ( 

if  (bndpts[curr_pt][2]  !=  bndpts[curr_pt  +  1][2])  { 
if  (bndpts[curr_pt][2j  =  bndpts[curr_pt  -79][2])  { 
regs(j][k]fO]  =  0.0; 
regs[j][k][l]  =  (float  )(curr_y  +  1.5); 
regs(j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

) 

else  if  (bndpts[curr_pt][2]  =  bndpts[curr_pt  -80][2])  { 
regs[j][k][0]  =  0.0; 

regs(j][k][l]  =  (float)(curr_y  +  0.5);  , 
regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

I 

else  { 

regsUmO]  =  0.0; 

regs[j][kj[l]  =  (float)(curr_y  -  0.5); 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 

k++; 


if  (hit_border)  | 
stop_tracing  =  TRUE; 

) 

else  { 

hit_border  =  TRUE; 
curr_pt  =  i; 

curr_x  =  bndpts|curr_pt][0]; 
curr_y  =  bndpts[curr_pt][l]; 
curr_dir  =  East; 


) 

else  if  (bndpts[curr_pt][2]  !=  bndpts[curr_pt  -  1][2])  { 

if  (bndpts[curr_pt][2]  ==  bndpts[curr_pt  -  81][2])  { 
regsUlfkHO]  =  0.0; 
regs[j][k][I]  =  (float )(curr_y  -  1.5); 
regsU]|k]l2]  =  (float)(bndptslcurr_pt]l2]); 
k++; 

1 

else  if  (bndpts[curr_pt][2]  ==  bndpts[curr_pt  -  80][2J)  { 
regs[j][k][0]  =  0.0; 
regs[j][k][l]  =  (float )(curr_y  -  0.5); 
regs[j][k][2]  =  (float  )(bndpts[curr_pt][2]); 
k++; 

I 

else  { 
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regs[j][k][OJ  =  0.0; 

regs[j][k][l]  =  (float  )(curr_y  +  0.5); 

regs{jJ[k](2J  =  (float  )(bndptsfcurr_pt][2]); 

k++; 

}  f*  3rd  level  else  */ 
if  (hit_border) 

stop_tracing  =  TRUE; 
else  | 

hit_border  =  TRUE; 
trace_counter  =  TRUE; 
curr_dir  =  South; 
curr_pt  =  i; 

curr_x  =  bndpts[curr_pt][0]; 
curr_y  =  bndpts[curr_pt][l]; 

} 

}  /*  2nd  level  else  */ 

else  if  (trace_counter)  { 

trace_border  =  TRUE; 
right_neighbor(); 
add_point(); 
curr_dir  =  North; 
move_curr_dir(); 

) 

else  if  (!trace_counter)  ( 
trace_border  =  TRUE; 
left_neighbor(); 
add_point(); 
curr_dir  =  South; 
tnove_curr_dir(); 

) 

)  /*  1st  level  if  */ 

if  (bndpts[curr_pt][0]  ==78)  ( 

if  (bndpts[curr__ptj[2]  !=  bndpts[curr_pt  +  1][2])  { 
if  (bndpts[curr_pt][2]  ==  bndpts[curr_pt  +  8 1 J [2] )  { 
regs[j][k]t0]  =  79.0; 
regs[j][k][l]  =  (float  )(curr_y  +  1.5); 
regs[j]lk][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

) 

else  if  (bndpts[cuir_pt][2]  =  bndpts[curr_pt  +  80][2J)  j 
regs[jJlk][0]  =  79.0; 
regs[j][kj[l]  =  (float  )(curr_y  +  0.5); 
regs[j][k][2]  =  (float  )(bndpts[curr_pt][2]); 
k++; 
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} 

else  { 

regs|j][k][0]  =  79.0; 
regs{j][k][l]  =  (float)(curr_y  -  0.5); 
regs(j][k][2]  =  (float)(bndpts[curr_ptl[2]); 
k++; 

) 

if  (hit_border) 

stop_tracing  =  TRUE; 
else  j 

hit_border  =  TRUE; 
trace_counter  =  TRUE; 
curr_dir  =  South; 
curr_pt  =  i; 

currx  =  bndpts[curr_pt][0]; 
curry  =  bndpts[cuir_pt][l]; 


else  if  (bndpts[curr_pt][2]  !=  bndpts[curr_pt  -  1][2J)  ( 
if  (bndpts[curr_pt][2]  =  bndpts[curr_pt  +  79][2J)  { 
regs[j][k][0]  =  79.0; 
regs[j]fk][l]  =  (float )(curr_y  -  1.5); 
regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

} 

else  if  (bndpts[curr_pt][2j  ==  bndpts[curr_pt  +  80][2J)  { 
regs(j][k][0]  =  79.0; 
regs[j][k][l]  =  (float )(curr_y  -  0.5); 
regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

) 

else  ( 

regs[j][k][0J  =  79.0; 
regs[jllk][l]  =  (float  )(curr_y  +  0.5); 
regs[j][k][2]  =  (float  )(bndpts[curr_pt][2]); 
k++; 

)  /*  3rd  level  else  */ 
if  (hit_border) 

stop_tracing  =  TRUE; 
else  { 

hit_border  =  TRUE; 
trace_counter  =  TRUE; 
curr_dir  =  South; 
curr_pt  =  i; 

curr_x  =  bndpts[curr_ptlfO]; 
curr_y  =  bndpts(curr_pt][  1  ]; 
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} 

}  /*  2nd  level  else  */ 

else  if  (trace_counter)  { 
trace_border  =  TRUE; 
right_neighbor(); 
add_point(); 
curr_dir  =  South; 
move_curr_dir(); 

} 

else  if  (!trace_counter)  ( 
trace_border  =  TRUE; 
left_neighbor(); 
add_point(); 
curr_dir  =  North; 
mo ve_curr_dir( ); 

) 

}  /*  1st  level  if  */ 

if  (bndpts[curr_pt][l]  —  1)  { 

if  (bndpts[curr_pt][2]  !=  bndpts[curr_pt  -80][2])  { 
if  (bndpts[curr_pt][2J  ==  bndpts[curr_pt  -  81][2])  ( 
regs[j][k][0]  =  (float )(curr_x  -  1.5); 
regs[j][k][l]  =  0.0; 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

I 

else  if  (bndpts[curr_pt][2]  ==  bndpts[curr_pt  -  1][2])  { 
regs|jj[k][0]  =  (float)(curr_x  -  0.5); 
regs[j]|k][l]  =  0.0; 

regs[j][k][2)  =  (float  )(bndpts[curr_pt][2]); 
k++; 

) 

else  ( 

regs[j][k][0]  =  (float  )(curr_x  +  0.5); 
regsljJlkJUl  =  0.0; 

regs[jj[k]12]  =  (float  )(bndpts[curr_pt][2]); 
k++; 

) 

if  (hit_border) 

stop_t  racing  =  TRUE; 
else  ( 

hit_border  =  TRUE; 
trace_counter  =  TRUE; 
curr_dir  =  South; 
curr_pt  =  i; 

curr_x  =  bndpts[curr_pt]f0]; 
curr_y  =  bndpts[curr_pt][l]; 
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} 

} 

else  if  (bndpts[curr_pt][2]  !=  bndpts[curr_pt  +  80][2])  { 
if  (bndpts[curr_pt][2]  =  bndpts[cun_pt  +  79][2])  { 
regs(JJ[k]lO]  =  (float  )(curr_x  +  1.5); 

regsUllkltl]  =  0.0; 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

} 

else  if  (bndpts[curr_pt][2]  =  bndpts[curr_pt  -  1][2])  { 
regs(j][k][0]  =  (float)(curr_x  +  0.5); 
regsUl[k)[l]  =  0.0; 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

} 

else  { 

regs[j][k][0]  =  (float)(curr_x  -  0.5); 
regs[j][k][l]  =  0.0; 

fegs[j][k][2]  =  (float  )(bndpts[curr_pt][2]); 
k++; 

)  /*  3rd  level  else  */ 
if  (hit_border)  { 

stop_tracing  =  TRUE; 

) 

else  ( 

hit_border  =  TRUE; 
trace_counter  =  TRUE; 
cun_dir  =  South; 
curr_pt  =  i; 

curr_x  =  bndpts[curr_pt][0]; 
curr_y  =  bndpts[curr_pt][l]; 

) 

}  /*  2nd  level  else  */ 

else  if  (trace_counter)  ( 
trace_border  =  TRUE; 
right_neighbor(); 
add_point( ); 
curr_dir  =  West; 
move_curr_dir( ); 

) 

else  if  (!trace_counter)  { 
trace_border  =  TRUE; 
left_neighbor(); 
add_point( ); 
curr_dir  =  East; 
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mo  ve_curr_dir( ); 

I 

}  f*  1st  level  else  */ 
if  (bndpts[curr_pt][l]  =  78)  { 

if  (bndpts[curr_pt][2]  !=  bndpts[curr_pt  -  80][2])  { 
if  (bndpts[curr_pt][2J  =  bndpts[curr_pt  -  79](2])  { 
regs[j][k]  [0]  =  (float  )(curr_x  -  1.5); 
regslj][k][l]  =  79.0; 

regs[j][kj[2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

I 

else  if  (bndpts[curr_pt][2]  =  bndpts[curr_pt  +  1][2])  { 
regs|j][k)[0]  =  (float)(curr_x  -  0.5); 
regs[j][k)[l]  =  79.0; 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

} 

else  ( 

regs[j][k][0]  =  (float  )(curr_x  +  0.5); 
regs[j][k][l]  =  79.0; 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 


if  (hit_border) 

stop_tracing  =  TRUE; 
else  ( 

hit_border  =  TRUE; 
trace_counter  =  TRUE; 
curr_dir  =  South; 
curr_pt  =  i: 

curr_x  =  bndptsfcurr_pt]fO]; 
curr_y  =  bndpts[curr_pt][l]; 


) 

else  if  (bndpts[curr_pt][2]  !=  bndpts[curr_pt  +  80][2])  { 
if  (bndpts[curr_pt][2]  ==  bndpts[curr_pt  +  81][2J)  { 
regs[jl[k]rO]  =  (float  )(curr_x  +  1.5); 
regs[j]fk][l]  =  79.0; 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

} 

else  if  (bndpts[curr_pt][2]  ==  bndpts[curr_pt  +  1  ][2])  { 
regs[j](k][0}  =  (fIoat)(curr_x  +  0.5); 
regs(j][k][  1  ]  =  79.0; 


52 


regs|j]M[2]  =  (float)(bndpts[cuir_pt][2J); 
k++; 

} 

else  { 

regs[j][k][0]  =  (float)(curr_x  -  0.5); 
regsUMl]  =  79.0; 

regs[j][k][2]  =  (float)(bndpts[curr_pt][2]); 
k++; 

)  /*  3rd  level  else  */ 

if  (hit_border) 

stop_tracing  =  TRUE; 
else  { 

hit_border  =  TRUE; 
trace_counter  =  TRUE; 
curr_dir  =  South; 
curr_pt  =  i; 

curr_x  =  bndpts[curr_pt][0]; 
curr_y  =  bndpts[curr_pt][l]; 

) 

J  /*  2nd  level  else  */ 
else  if  (trace_counter)  { 
trace_border  =  TRUE; 
right_neighbor(); 
add_point(); 
curr_dir  =  East; 
move_curr_dirO; 

) 

else  if  (!trace_counter)  ( 
trace_border  =  TRUE; 
left_neighbor(); 
add_pointO; 
curr_dir  =  West; 
move_curr_dir( ); 


|  /*  1st  level  if  */ 
hit_border  =  TRUE; 


}  /*  function  connect_border  */ 


/«.**t**  +  +  ;f*  +  t****1<*  +  **  +  **  +  *******^  +  ***  +  *  +  +  +  *4,*  +  **  +  *  +  +  <,***  +  4[^4,t  +  +  **4i*  +  *  +  *  ++ J 


I*  Inserts  a  point  to  the  current  region  output  database,  and  increments  */ 
/*  the  point  count.  */ 
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add _point() 

{ 

regsfj]fk][0]  =  (float)  (curr_x  +  neighbor_x)/2; 
regs(j][k][l]  =  (float)  (curr__y  +  neighbor_y)/2; 
regsjj)[kj[2]  =  (float)  (bndpts[curr_pt][2]); 
k++;  /*  increment  next  point  to  be  added  to  the  region  */ 


}  /*  function  add_point  */ 


f*  *********************************************************************4 <*/ 


/*  Inserts  a  point  to  the  current  region  output  database,  and  increments  */ 
/*  the  point  count.  This  is  used  when  the  region  is  closed,  and  the  */ 
/*  first  point  is  the  same  as  the  last  point.  *1 


duplicate_first_pt() 

I 

regs[j]P0r0]  =  regs(j3f0][0]; 
regs|j][k][l]  =  regs(j][0][l]; 

k++;  /*  increment  next  point  to  be  added  to  the  region  */ 

)  /*  function  duplicate_first_pt  */ 

^******************************************4 '*****************************/ 


/*  Gets  the  x,y  of  the  point  to  the  left  of  the  curr  point.  */ 

/*  If  we  are  tracing  counterclockwise, :we  get  the  x,y  of  point  to  the  right  */ 

left_neighbor( ) 


switch  (curr_dir) 

( 

case  North  : 


case  NEast  : 


case  East  : 


case  SEast  : 


case  South  : 


neighbor.^  =  curr_x  -  1; 
neighbor_y  =  curr_y; 
break; 

neighbor_x  =  curr_x  -1; 
neighbor_y  =  curr_y  +1; 
break; 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  +  1; 
break; 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  +  1; 
break; 

neighbor_x  =  curr_x  +1; 


54 


neighbor_y  =  curr_y; 
break; 


case  SWest 

neighbors  =  curr_x  +  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case  West  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  -  1; 
break; 

case  NWest 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  -  1; 
break; 

default: 

printf("Error  curr-dir  !=  0-7"); 
break; 

}  /*  switch 

*/ 

}  /*  function  */ 

/*  Gets  the  x,y  of  the  left  rear  neigbor  of  the  current  point.  */ 
left_rear_neighbor() 


switch  (curr  dir) 

I 

case  North  : 

neighbor_x  =  curr_x  -  1; 
neighbor _y  =  curr_y  -  1 ; 
break; 

case  NEast  : 

neighbor_x  =  curr_x  -1; 
neighbor_y  =  curr_y; 
break; 

case  East  : 

neighbor_x  =  curr_x  -  1; 
neighbor_v  =  curr_y  +  1; 
break; 

case  SEast  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  +  1; 
break; 

case  South  : 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  +  1; 
break; 

case  SWest  : 

neighbor_x  =  curr_x  +  1 
neighbor_y  =  curr_y; 
break; 

case  West  : 

neighbor_x  =  curr_x  +  1 
neighbor_y  =  curr_y  -  1; 
break; 
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case  NWest  :  neighbor_x  =  curr_x; 

neighbor_y  =  curr_y  -  1; 
break; 

default:  printf("Error  curr-dir  !=  0-7"); 

break; 

}  /*  switch  */ 


)  /*  function  /* 


^* *****************************************************************  j 

/*  Gets  the  x,y  of  the  left  front  neighbor  of  the  current  point.  */ 


left_fr  ont_neighbor( ) 

{ 

switch  (curr_dir) 


case  North  : 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  +  1; 
break; 

case  NEast  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  +1; 
break; 

case  East  : 

neir'ibor_x  =  curr_x  +  1; 
neighborly  =  curr_y  +  1; 
break; 

case  SEast  : 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y; 
break: 

case  South  : 

neighbor_x  =  curr_x  +1; 
neighbor_y  =  curr_y  -  1; 
break; 

case  SWest  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  -  1; 
break; 

case  West  : 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case  NWest 

:  neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y; 
break; 

default: 

printfC'Error  curr-dir  !=  0-7"); 
break; 

)  /*  switch  *1 


)  /*  function  */ 
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I*  Gets  the  x,y  of  the  front  neighbor  of  the  current  point.  */ 
front_neighbor() 


if  (!  trace_counter)  ( 
switch  (curr_dir) 

{ 


case  North  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  +  1; 
break; 

case  NEast  : 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  +  1; 
break; 

case  East  ; 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y; 
break; 

case  SEast  : 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case  South  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  -  1; 
break; 

case  SWest  : 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case  West  : 

neighbor_x  =  curr_x  -  1: 
neighbor_y  =  curr_y; 
break; 

case  NWest  : 

neighbors  =  curr_x  -  1; 
neighbor_y  =  curr_y  +  1: 
break; 

default: 

)  /*  switch  */ 

)  /*  if  trace  left  */ 
else  { 

switch  (curr_dir) 

printf("Error  curr-dir  !=  0-7 
break; 

t 

case  North  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  +  1; 
break; 

case  NEast  : 

neighbor_x  =  curr^x  +  1; 
neighbor_y  =  curr_y  +  1; 
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break; 

case  East  :  neighbor_x  =  curr_x  +  1; 

neighborly  =  curr_y; 
break; 

case  SEast  :  neighbor_x  =  curr_x  +  1; 

neighbor_y  =  curr_y  -  1; 
break; 

case  South  :  neighbor_x  =  curr_x; 

neighbor_y  =  curr_y  -  1; 
break; 

case  SWest  :  neighbor_x  =  curr_x  -  1; 

neighborly  =  curr_y  -  1; 
break; 

case  West  :  neighbor_x  =  curr_x  -  1; 

neighbor_y  =  curr_y; 
break; 

case  NWest  :  neighbor_x  =  curr_x  -  1; 

neighbor_y  =  curr_y  +  1; 
break; 

default:  printf("Error  curr-dir  !=  0-7"); 

break; 

}  /♦  switch  */ 

J  /*  else  trace  right  */ 

)  /*  function  */ 

********  *#  ********  a***********  ****************  ******************** ******/ 

/*  Gets  the  x,y  of  the  point  to  the  left  of  the  curr  point.  */ 

/*  If  we  are  tracing  counterclockwise,  we  get  the  x,y  of  point  to  the  right  */ 

right_neighbor( ) 

I 

switch  (curr_dir) 

( 

case  North  :  neighbor_x  =  curr_x  +  1; 

neighbor_y  =  curr_y; 
break; 

case  NEast  :  neighbor_x  =  curr_x  +1; 

neighbor_y  =  curr_y  -1; 
break; 

case  East  :  neighbor_x  =  curr_x; 

neighbor_y  =  curr_y  -  1; 
break; 

case  SEast  :  neighbor_x  =  curr_x  -  1; 

neighbor_y  =  curr_y  -  1; 
break; 
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case  South  : 

neighbor,*  =  curr_x  -1; 
neighbor_y  =  curr_y; 
break; 

case  SWest  : 

neighbor_x  =  curr_x  -  1; 
neighborly  =  curr_y  +  1; 
break; 

case  West  : 

neigh  bor_x  =  curr_x; 
neighbor_y  =  curr_y  +  1; 
break; 

case  NWest 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  +  1; 
break; 

default: 

printf("Error  curr-dir  !=  0-7"); 

break; 

)  /*  switch  */ 

)  /*  function  */ 


p*  ************************************************************************ ^ 


/*  Gets  the  x,y  of  the  point  to  the  right  front  of  the  curr  point.  */ 
right_front_neighbor( ) 

I 

switch  (curr_dir) 


case 

North  : 

neighbor_x  =  curr_x  +  1: 
neighbor_y  =  curr_v  +  1; 
break; 

case 

NEast  : 

neighbor_x  =  curr_x  +  1 
neighbor_y  =  cun_y; 
break; 

case 

East  : 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case 

SEast  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  -  1; 
break; 

case 

South  : 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case 

SWest  : 

neighbor_x  =  curr_x  -  1 
neighbor_y  =  curr_v; 
break; 

case 

West  : 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  +  1; 
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break; 

case  NWest  :  neighbor_x  =  curr_x; 

neighbor_v  =  curr_y  +  1; 
break; 

default:  printf("Error  curr-dir  !=  0-7"); 

break; 

}  f*  switch  */ 

}  /*  function  */ 


f*  Gets  the  x,y  of  the  point  to  the  right  rear  of  the  curr  point.  */ 

right_rear_neighbor( ) 

I 

switch  (curr-dir) 


case 

North  : 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case 

NEast  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  -  1; 
break; 

case 

East  : 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  -  1; 
break; 

case 

SEast  : 

neighbor_x  =  curr_x  -  1; 
neighboi_y  =  curr_y; 
break; 

case 

South  : 

neighbor_x  =  curr_x  -  1; 
neighbor_y  =  curr_y  +  1; 
break; 

case 

SWest  : 

neighbor_x  =  curr_x; 
neighbor_y  =  curr_y  +  1; 
break; 

case 

West  : 

neighbor_x  =  curr_x  +  1; 
neighbor_y  =  curr_y  +  1; 
break; 

case 

NWest 

neighbor_x  =  curr_x  +  1 
neigh  bor_y  =  curr_y; 
break; 

default: 

printfC'Error  curr-dir  !=  0-’ 

break; 

)  /*  switch  */ 

)  /*  function  */ 
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/*  Marks  the  current  point  as  used  in  a  region,  and  changes  the  current  */ 
f*  point  in  the  direction  of  the  current  direction.  Changes  the  curr_x,  */ 

/*  and  curr_y  to  correspond  with  the  new  current  point.  */ 


move_curr_dir() 

I 

bndpts[curr_pt]f4]  =  TRUE; 
switch  (curr_dir) 

I 

case  North  :  curr_pt  =  curr_pt  +  1; 
break; 

case  NEast  :  curr_pt  =  curr_pt  +  81; 
break; 

case  East  :  curr_pt  =  curr_pt  +  80; 

break; 


case  SEast  :  curr_pt  =  curr_pt  +  79; 
break; 

case  South  :  currpt  =  curr_pt  -  1; 
break; 

case  SWest  :  curr_pt  =  curr_pt  -  81; 
break; 

case  West  :  curr  pt  =  curr_pt  -  80; 

break; 

case  NWest  :  curr_pt  =  curr_pt  -  79; 
break; 

default  ;  printf("Curr_dir  !=  0-7”); 
break; 

} 

curr_x  =  bndpts[curr_pt][0]; 

curr_y  =  bndpts[curr_pt][l]; 

}  /*  function  move_curr_dir  */ 

^*  *****%********★*******************************************************/ 


print_output_regions( ) 

I 

FILE  *fdout; 

FILE  *foutfig: 
unsigned  short  r; 
unsigned  short  s; 
int  temp; 
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fdout  =  fopen(outfile,  "w+"); 
foutfig  =  fopen(outfig,  ”w+"); 
if  (fdout  ==  NULL)  { 
printfO'Have  file  opening  problem"); 

} 

if  (foutfig  =  NULL)  { 
printfO'Have  file  opening  problem"); 

) 

/*  Writes  to  file  and  screen  the  x  and  y  value  for  each  of  the  regions.  */ 
f*  The  stopping  conditions  of  the  loops  are  when  the  number  of  regions  */ 
f*  generated  is  reached,  and  for  each  region  when  the  flag  9999  is  reached.  */ 
/*  This  function  assumes  it  knows  the  number  of  regions  (j  -1)  and  the  */ 
/*  array  element  after  the  last  point  in  each  region  is  flagged.  */ 


fprintf(fdout. "module  regions. W); 
f^rintf(  fdout, "/*$eject*An"); 
fj>rintf(  fdout, "body  .W ); 
for  (r  =  0;  r  <  j;  r++)  { 

fprintf(fdout,"reg(%hu,  ["/); 
s  =  0; 
do  { 

/*  Convert  to  inches  only  for  the  foutfig  file.  */ 

printf("%4hu  %4f  %4f\n",  r,regs[r]fslfO].regs[r][s][l]); 
fprintf(fdout,"[%.lfeO,  %.lfe0]'\regs[rj[s][0], 
regs[r][sj[lj); 

regsfr][s]fO]  =  (regsfr][s][0]  /  10); 
regs[r][s]fl]  =  (regs[r][s]f  1]  /  10); 
fprintf(foutfig,"%4f  %4fSn  'Vegs[r][s][0], 
regs[r][s][l|); 

fprintf(foutfig.”%4f  %4f\n  ",regs[r][s][0], 
regs[r][s]ll]); 
s++; 

if  (regs[r][s][0]  !=  9999.000000)  f 
fprintf(fdout,",\n" ); 

) 

)  while  ((regs[r][s][01  !=  9999.000000)  II 
(regs[r][s][l]  !=  9999.0000000)); 
temp  =  (int)(regsfr](01f2]); 
fprintf(fdout,"J,%hu,%hu).Vi",temp,s); 

)  /*  for  r  */ 

fprintf(  fdout .  "numreg  ions(  %hu )  ,Nn  "  ,r ) ; 

f))rintf(fdout."endmod.\n"); 

fclose(fdout); 
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fclose(foutfig); 

}  /*  Print  the  output  of  the  regions.  */ 

j*  F.tvt  ofBNDUTIL.C  *****************************************************•/ 
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File  header.h 

/*  This  file  sets  constants  and  determines  the  path  for  the  input  and  */ 
/*  output  file  for  the  bndpts.c  program.  */ 

#define  infile  "/work/wade/Thesis/C/rawdata2.h" 

#defme  outftle  "/work/wade/Thesis/C/regiondata2.pro" 

#define  outfig  "/work/wade/Thesis/C/regiondata2.fig" 

#define  boolean  int 
#define  FALSE  0 
#define  TRUE  ! (FALSE) 

#define  North  0 
#define  NEast  1 
#define  East  2 
#define  SEast  3 
#define  South  4 
#define  SWest  5 
#define  West  6 
#define  NWest  7 

/********END  of  file  header.h  ***************************/ 
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File  vdb2.c 

/*  This  file  is  the  main  program  for  getting  a  1  Km  by  1  Km  grid  square  */ 
l*  from  the  DTED  data  base.  The  program  was  run  on  the  Silicon  Graphics  */ 
t*  Computer.  The  program  requires  file  "files.h"  to  run.  The  program  is  */ 

I*  a  variant  of  a  widely  used  program  in  the  department  to  get  the  data  */ 

#include  "stdio.h" 

#include  "ctype.h" 

#include  "math.h" 

#include  "files.h" 

char  infile[50]=  MASTER_DMA_DTED_FI LE ; 
int  fdin; 

main() 

I 

int  fdout,x.z.r,c,i,j,utmx,utmz; 

int  off,ewerror=l,nserror=l,doswap=l; 

char  s,swap,outfile[50],utmew|5j,utmns[5],edbase[20402],temp; 

unsigned  short  raw .elev.veg, massaged; 

shon  MAXUTMEW=659,MAXUTMNS=849,DATAPTS=1 1; 

systemf'clear"); 

printffTHIS  PROGRAM  WILL  CREATE  A  TERRAIN  DATABASEAn”); 
printf("THE  SOURCE  DATABASE  IS  A  DEFENSE  MAPPING  AGENCY 
DIGITALVi"); 

printf( "TERRAIN  ELEVATION  FILE  FOR  A  36  KM  BY  35  KM  REGION  OF\n"); 

print f(" FORT  HUNTER  LIGGET,  CA.  AND  VICINITY  An"); 

pnntfC'THE  OUTPUT  FILE  IS  A  1  KM  BY  1  KM  SUBSET  OF  THE  ENTER  EAn"); 

printffREGION.  AND  WILL  BE  STORED  IN  THE  FORMAT  REQUIREDNn"); 
printf("  100.0  meter  resolution\n”). 

printft,"  80  x  80  data  pointsNn"); 

print f("  storage  in  z  major  orderVi"); 

pnntfC'YOU  MUST  ENTER  THE  LNM  COORDINATES  OF  THE  SOUTHWESTVi"); 

printfCCORNER  OF  THE  SUBSET  REGION  YOU  WANT  EXTRACTED  An"); 
printfC  ENTER  A  CARRIAGE  RETURN  TO  CONTINUE...”); 


while(ewerror) 
systemt  "clear”); 
printft" 
printf(" 
prmtf(" 


VALID  EW  COORDINATES  ARE  IN  THE  RANGES  An"); 
EAST- WEST  (UTM  EW):  410  to  %cfnV\MAXUTMEW); 
*******************v-); 
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printfC 

printfC 

printfC 

printfC 

printfC 

printf(" 

printfC 

printfC 

printfC 

printfC 

printfC 

printf(" 


*  **T); 

*  *\n"); 

*  N  *\n"); 

*  I  *\n"); 

*  W - E  *\n"); 

*  I  *«"); 

*  S  *\n"); 

*  *Sn"); 

*  *\n"); 
*******************^"^ 

I  NO; 

410  %dSn",MAXUTMEW); 


printf("\n\nENTER  the  UTM  EW  coordinate  of  the  southwest  comeiNn"); 
printfC  (  Enter  an  integer  value  X:  410  <=  X  <=  %d  )Nn"fMAXUTMEW); 
printfC  X  ?  =>"); 

scanf(  "%d"  ,&utmx ) ; 

if(  (utmx>=410)  &&  (utmx<=MAXUTMEW)  )  ewerror=0; 


while(nserror)  { 
system(  "clear"); 

printf("  VALID  NS  COORDINATES  ARE  IN  THE  RANGES:V'); 

printfC  NORTH-SOUTH  (UTM  NS):  600  to  %cLnW\MAXUTMNS); 

printfC  %3d-  *******************W\MAXUTMNS); 

printfC  *  *W); 

printfC  *  *V'); 

printf("  *  N  *\n”); 

printfC  *  I  *W'); 

printfC  *  W- — E  *Vn"); 

printfC  *  I  *\n"); 

printf("  *  S  *\n"); 

printf("  *  *\n"); 

printfC  *  *\n"); 

printf("  600-  fr******************^")- 

printf(”\n\nENTER  the  UTM  NS  coordinate  of  the  southwest  comerVi"); 

printfC  (  Enter  an  integer  value  Y:  600  <=  Z  <=  %d  )Nn">lAXUTMNS); 

printfC  Z  ?  ==>"); 

scanf("%d",&utmz); 

if(  (utmz>=600)&&(utmz<=MAXUTMNS)  )  nserror=0; 


systemCclear"): 

printfC  DO  YOU  WISH  ELEVATION  DATA  WORDS  TO  HAVE  THEIR  BYTES 
SWAPPEDTNn" ); 

printfC  ENTER  y  (  YES:  SWAP  BYTES  )V); 

printfC  ENTER  ’n’  (  NO:  NO  BYTE  SWAP  )W>: 

printfC  SWAP  BYTES  ==>  ?"); 
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swap  =  toupper(getcharO); 

while(  (swap  !=  ’Y’)  &&  (swap  !=  ’N’)  )  swap=touppcr(getchar()); 

if(swap=’N’)  doswap=0; 

system("clear"); 

sprintf (utmew , "  %  3  d "  ,utmx ) ; 

sprintf(utmns,"%3d",utmz); 

strcpy(outfile,OUTPUTFILE); 

if(doswap)  strcat(outfile,"swap"); 

printf("\nCreating  database  for  a  1km  x  1km  region,  southwest  comer"); 
printf("Nnat  %d  -  %d.  Database  will  consist  of  elevation  data’’,utmx,utmz); 
printfC'Nnfor  %d  x  %d  points  at  12.5  meter  resolution.", DATAPTS.DAT APTS); 
printf("ViDatabase  filename  is  %s\n",outfile); 
fdin  =  open(infile,0); 
fdout=  creat(outfile,0644); 


r  =  (utmz-600)*8; 
c  =  (utmx-410)*8; 

Iseek(fdin,offset(r,c),0); 
for  (i  =  0;  i<  80;  i++){ 
for(j  =  0;  j<  80;  j++)| 

read(fdin,&raw,2); 

veg  =  ((unsigned  short)(raw  &  OxeOOO)  »  13>+48; 

write(fdout,&veg,2); 

printf("%c",veg); 

) 

printfCVi"); 

) 

close(fdout); 

close(fdin): 

|  /*  main  */ 

/*  This  calculates  the  startpoint  for  the  gridsquare  within  the  */ 

/*  36  Km  by  35  Km  database.The  6400  represents  the  number  of  data  points  */ 
/*  per  grid  square.  The  35  is  the  number  of  grid  squares  in  the  north/south  */ 

/*  dnection.  */ 

offset(r,c) 
int  r,c; 

( 

return  (  2  *  ( 

(6400  *  35  *  (int  )(c/80)) 

+  (6400  *  (int )( r/80 ) ) 

+  (  80  *  (c%80)) 

+  (  (r%80)) 
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) 


); 

} 

/+***+** gn(j  gjg  VDB2.C  ***’*,‘*'’*'’*'’*''^’^^**1*,,*'*'*'^,*,’*,*J*'**1*'****'*'’*‘’*'*^***'*'^*/ 

File  files. h 

/*  This  file  sets  the  paths  required  for  input/output  to  vdb2.c  */ 

#define  MASTER_DMA_DTED_FILE  "/usr/work/cdec/DTED/terrain.dat" 
#define  OUTPUTF1LE  "usr/work/wade/thesis/rawdata.h" 

#define  PRINTFILE  "usr/work/wade/thesis/rawdata.p" 

^♦******£0^ j  q£  files. h  ’fr**************************************/ 
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APPENDIX  C  SOURCE  CODE  FOR  THE  SPLIT-AND-MERGE  PHASE 


module  split -and -merge. 

impoit(reg  /4). 

/*$eject*/ 

body. 

dynamic(min_index/l ). 
dynamic(max_index/l ). 
dynamic(segment/4 ) . 
dynamic(newregions/3 ). 
dynamic(outputlist/l ). 
dynamic(intermedlist/l ). 
dynamic(currentregion/l ). 
dynamic(flag/l ). 
dynamic(pairlistin/l ). 
dynamic(pairlistout  1/1 ). 
dynamic(pairlistout2/l ). 
dynamic(sproditem/0). 
dynamic(ssumitem/0). 
dynamic(incr _global/2). 
dynamic(vprodout/l ). 
dynamic(vprodinl/l ). 
dynamic! vprodin2/l ). 
dynamic(sumuplist/l ). 
dynamic(sumupsum/l ). 

split_threshold(  100e-2) 
mergethresholdl  100e-2). 

gol 

set_state(global_stack,30000), 

set_state(main_stack,  1 0000), 

system(compress_stacks), 

system(garbage_collection), 

display_statistics, 

handle_adj_regs(reg,newregions), 

split_merge(  reg.newregions), 

print_output. 

nl. 
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display_statistics 
nl,  stars, 

write_tab(18),  write("SYSTHM  STATUS”),nl, 
state(cpu_time  ,X) , 

write("cpu  time  =  "),  write(X),write("  msec"),  nl, 
state(main_stack,  [U,C]), 
write("main  stack  used  =  "),  write(U),nI, 
state(global_stack,  [UG,CG]), 
write("global  stack  used  =  "),  write(UG),nl, 
state(statement_table,  IU 1  ,C  1  ] ), 
write( "statement  table  used  =  "),  write(Ul),  nl, 
stars,  nl. 

stars  -  writeC***********************************************************") 
nl. 

print_output 

newregions(Regnumber,NewPointList, Value), 

write("The  "),write(Regnumber),write("  Region  is:  "),nl, 

write(NewPointList),nl, 

write("The  Value  is:  "),  write(Value),nl,nl, 

fail. 

print_output. 

handle_adj_regs(INNAME,  OUTNAME)  :- 

del_all_statements(OUTNAME/3  ),handle_adj  1  (INNAME,OUTN  AME), ! . 

handle_adjl  (INNAME, OUTNAME)  :- 

get_a_regionfINNAME.Regl  ,PL1  .VI  .NP1 ), 
get_a_region(INNAME,Reg2,PL2,V2,NP2), 

Regl  =/=  Reg2, 
display_statistics. 
check_if_adj(PLl.PL2,AdjPts),!, 
display_statistics. 

handle_adj2(OLTN  AME, Reg  1  ,PL  1 ,  V 1  ,NP  1  ,Reg2,PL2,V2,NP2,AdjPts ). 
fail. 

handle_adj  1  (INN  AME, OUTNAME;. 

check_if_adj(PLl.PL2.AdjPts)  :- 
real_uitersection(PLl,PL2,AdjPts), 
length(AdjPts.LA),  write(LA),nl,LA  >  0, 


handle_adj2(OUTNAME,Regl,PLl,Vl,NPl.Reg2,PL2,V2.NP2.AdjPts)  :- 
display_statistics, 

handle_adj3(Regl  ,PL1  ,NP1  .AdjPts.NewPLl ), 
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OP  =..  [OUTNAME,  Regl,  NewPLl,  VI], 
asseitz(OP ), 

handle_adj3(Reg2,PL2,V2,NP2,AdjPts,NewPL2), 

OQ  [OUTNAME,  Reg2,  NewPL2,  V2], 

assertz(OQ), 

fail. 

handle_adj2(INNAME, OUTNAME)  !. 

handle_adj3(Reg.PL,NPAdjPts,NewPL) 
write("inside  adj3"),nl, 

set_global(currentregion,PL), 
sct_global(outputlist,[]), 
set _global(intermedlist,[]), 
set_global(flag,l), 

find_sublist_indices(PL,N  1  ,N2,AdjPts), 

get_sublist(PL,l  ,N1  ,FPL1  ),get_sublist(PL,N2,NPJBPLl ), 

length(AdjPts,LA), 

asserta(segment(Reg,Nl,N2,100)). 

split_into_segments, 

asserta(segment(Reg,  1  .N 1 , 1 00) ), 

asserta(segment(Reg,N2,NP,100)), 

spl  it_into_segments, 

set_global(flag.l).set_global(min_index,l  ),set  _global(max_index,NP), 

adj_merge(Reg.N  1  ,N2), 

build_intenned_list(Reg), 

intemiedlist(NewlndexList),  write(NewIndexList),nl, 
length(NewlndexList,NewNum),  write("new  list  length  is  "), 
write(NewNum),  nl.  buildnewlist(NewIndexList.NewPL), 
wr:te("The  Final  Point  List  is  "),  write(NewPL),  nl, 
display_segments_asserted,! . 

adj_merge(R.Pl,P2) 

display_segments_  asserted, nl. 
flag(l),  set _global(flag.O), 
doall(adj_merge_segment(R,P  1  ,P2 ) ). 
adj_merge(R,Pl,P2). 
adj_merge(R,Pl  ,P2 ). 

adj_merge_segment(R,Pl,P2) 

segment(R.A.B,FAB),segment(R,C,D,FCD).  B  =  C,  B  =/=  PI. 

B  =/=  P2,  not  segment(RtA,D,FAD). 
mcrge_segment2(R,A,B,C,D,FAB,FCD). 

adj_merge_segment(R.Pl,P2) 

segment(R.A.B,FAB  ).segment(R.C,D,FCD),  min_  index  (MIN), 
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max_index(MAX),C  =  MIN,  B  =  MAX,  B  =/=  PI,  C  =/=  PI, 

B  =/=  P2,  C  =/=  P2, 
merge_segment3(R,A,B,C,D,FAB,FCD). 

split_merge(INNAME,  OUTNAME) 
sp!it_merge2(INNAME,0UTNAME). 

sp!it_merge2(INNAME,  OUTNAME) 

get_a_region(INN  AME  .RegNumber  .PointList  .Value  .NumPoints) , 
split_merge3(RegNumber,PointList,NumPoints), 
set_gIobaI(flag,l ),  buiId_intermed_list(RegNumber), 
intermedlist(NewIndexList),  write(NewIndexList),nl, 

length(NewIndexList,NewNum),  write("new  list  length  is  "),write(NewNum), 

buildnewlist(NewIndexList,NewPointList), 

write("The  Final  Point  List  is  "),write(NewPointList),nl, 

OP  =..  [OUTNAME, RegN umber, NewPointList, Value], 

asseitz(OP), 

fail. 

split_merge2(INNAME,  OUTNAME)  !. 

get_a_region(INNAME,RegNumber.PointList, Value, NumPoints) 
not  segment(RegNumber,PointList, Value, NumPoints), 

Q  =..  [INNAME,RegNumber,PointList,Value,NuniPoints),  call(Q). 


/*  This  is  the  main  control  structure  for  the  split  and  merge  algorithm.  */ 

split_merge3(R,PointList.NumPoints) 
set  __global(currentregion,PointList), 
set_gl  obal  ( outpu  tl  ist .  [  ] ) , 
set  _global(intermedlist,[]), 
set _global(flag.l ), 

asserta(segment(R,l, NumPoints,  100)). 
split_into_segments, 
display_statistics, 
set  _global(flag,l), 
set_global(min_index,l ), 
set_global(max_indexNumPoints), 
merge_segments(R ), 
display_segments_asserted, ! . 

split_into_segments 

display_segments_assened,nl. 

flag(l).  set_global(flag.O).  doalI(spIit_into_segment), 

split_into_segments. 
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split_into_segments. 


split_into_segment 

segment(R,N  1  ,N2  ,F)  ,split_into_segment  2(R,N  1  ,N2,F). 

split_into_segment2(R,Nl,N2fF) 

split_threshold(T), 

F  >  T,  average(Nl,N2,N3),  retract(segment(R,Nl>N2,F)),  !, 
display_statistics, 

writefThe  first  sublist  indices  and  fit  are:  "), 
write(N  1 ), write("  "),write(N3),write("  "),nl, 
find_fit(R.N  1  ,N3  ,F  1 3 ),  asserta(segment(R,N  1  ,N3 , FI  3 )), 
write(F13).nl, 
display  _statistics, 

write("The  second  sublist  indices  and  fit  are:  "), 

write(N3), write!”  "),write(N2),write("  "), 

find_fit(R,N3,N2,F32),asserta(segment(R,N3,N2,F32)), 

write(F32),nl. 

set  _global(flag.l),!. 

merge_segments(R) 

display_segments„asserted,nl, 
flag(l),  set _global(flag.0), 
doall(merge_segment(R ) ) , 
merge_segments(R ) . 
merge_segments(R ). 

merge_segment(R) 

segment(R,Nl,N2,F12),segment(R,N3,N4,F34),  N2  =  N3, 
not  segment(R,Nl,N4,F3), 

nl, write!"  Nl=  ”),write(Nl ),write("  N2  =  "),write(N2), 
write("  N3=  "),write(N3).write("  N4  =  "),write(N4),nl, 
merge_segment2(R.N  1  ,N2.N3,N4,F1 2.F34). 
merge_segment(R) 

segment(R,Nl  .N2,F12),segment(R,N3,N4.F34), 
rmn_  index  (MIN  ],max_index(MAX ), 

N3  =  MIN, 

N2  =  MAX. 

write!"  Nl=  "),write(Nl),write("  N2  =  ”),write(N2), 
write!"  N3=  ").write(N3), write!"  N4  =  "),write!N4),nl, 
merge_segment3(R.N  1  ,N2,N3,N4,F  1 2,F34). 


merge_segment2(R.Nl.N2.N3,N4.F12.F34) 

display_statistics. 


73 


find_fit(R,Nl,N4,F14), 
writc("Fit  is  "),write(F14),nl, 
mergethreshold(T),  FI  4  <=  T, 
asserta(scgment(R,Nl  ,N4,F14)), 
retract(segment(R,N  1  ,N2,F12)), 
retract(segment(R,N3tN4,F34)), 
sct_global(flag,l),  !. 

merge_segment3(R,Nl,N2,N3,N4,F12,F34) 
display_statistics, 
find_wrap_fit(R,N  1  ,N4,F14), 
write("Fit  is  :  "),write(F14),nl, 
mergethreshold(T),  F14  <=  T, 
asserta(segment(R,N  1  ,N4,F  1 4)), 
retract(segment(R,N  1  ,N2,  F 12)), 
retract(segment(R,N3,N4,F34)), 
set_global(min_index,N4), 
set_global(max_index  ,N  1 ), 
set_global(flag,I),  !. 


display_segments_asserted 

segment(A,B,C,D),  write("Reg=  "),  write(A), 

write("  Indices  =  "),  write(B), 

write("  "),  write(C), 

write("  Fit  =  ").  write(D), 

nl,  fail. 

display_segments_asserted  nl. 

average(Nl.N2.N3) 

NP1  is  Nl  +  1,  NP1  =/=  N2, 

N3  is  ((N2  -  Nl)  div  2)  +  Nl,!. 

find_fit(R,A,B,Fit) 

currentregion(X), 

A  <  B, 

get_sublist(X.A.B.SubList),nl, 
endptlsline(SubList,[Cl,C2,C3]), 
lslinefit(SubList,[C  1  ,C2,C3]  ,Fit ). 

find_wrap_fit(R.A,B,Fit) 

currentregion(X), 

A  >  B, 

get_rest_list(X,A,SubListl ), 
get_sublist(X.l  ,B,SubList2), 
append(SubListl.SubList2,SubList). 
endptlsline(SubList,lC  1  ,C2,C3] ). 
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lslinefit(Sub  List,  [C 1  ,C2,C3] ,  Fit). 


build_inteimed_list(R) 
flag(  1  ),set_global(flag,0), 
doail(gct_segments(R)), 
intermedlist(IL), 

unduplicate(IL,NewIL),  sort(NewIL,NNewIL), 
set_global(intcrmedlist  ,NNe  wIL). 

build_intermed_Iist(R)  !. 

get_segments(R) 

segment(R,N  1  ,N2,Fit), 
get_segments2(R,N  1  ,N2,Fit). 

get_segmcnts2(R,N  1  ,N2JFit):- 
cons_global(intermedlist,Nl ), 
cons_global(intermedlist,N2), 
set_global(flag.l),!. 


buildnewlist(NewIndexList,NewPointList) 
currentregion(PointList ), 
set_global  (output!  ist ,  [  J ) , 
buildnewlist2(PointList,NewIndexList), 

outputlist(RNewPointList),  reverse(RNewPointList,NewPointList),! . 

buildnewlist2(PointList,G ). 
bui]dne\vlist2^PointList.rN!BNewIndexList]) 
get_point(Point  List  ,N  .Point ), 
cons_global(outputlist.Point), 
buiidnewlist2(PointList.BNewIndexList). 


get_point([AIBPointList],N.A)  :* 

N  =  1,  !. 

gct_point([AIBPointList] J^.Point) 

NC  is  N  -  1, 

get_point(BPointList,NC, Point),  !. 

/*  get_sublist  is  called  when  List,  Nl,  and  N2  are  bound:  yields  SubList  */ 
get_sublist(List.Nl,N2,SubList) 
length(List,LL),  LL  >=  2, 
get_rest_list(List,Nl,Ll ),(NC2  is  N2  -  Nl  +  2), 
get_rest_list(Ll.NC2  ,L2), 
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append(SubList,L2,Ll ). 


get_rest_l ist(Re  st  List  ,N,RestList ) 

N  =  1,  !. 

get_rest_list([AIL],N,Restlist)  > 

NC  is  N  -1, 

gct_rcst_list(L,NCrRestList),!. 

/*  flnd_sublist_indices  is  called  when  List  and  SubList  are  bound.  */ 

find_sublist_indices(List,Nl  ,N2,[FirstlB  SubList]) 
last(B  SubList  ,Last ) , 
find_sublist_indices2(List,l,N3, First), 
find_sublist_indices2(  List ,  1  ,N4,Last ) , 
order_indices(N3,N4,N  1  ,N2),  !. 

find_sublist_indices2([[Xl  ,Y13IL],NC,NC,[X2,Y2]):-  closer(Xl  ,X2),  closer(Y2,Y2), 

find_subIist_indices2([AIL],NC,N,Point) 

NC2  is  NC  +  l,find_sublist_indices2(L,NC2,N,Point),  !. 

order_indices(Nl,N2,Nl,N2)  N1  <=  N2,!. 
order_indices(Nl,N2,N2,Nl)  !. 

*doall(P)  not  alltried(P). 

*alltried(P)  call(P),  fail. 
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module  linear  least-squares. 

*split_pairs(  L ,L  1  ,L2 )  set_global(paiiiistin,L),  set _global(pairlistoutl,0), 

set_global(pairlistout2,Q),  iteratelist(split_pair(L),pairlistin), 
pairlistoutl(RLl),  pairiistout2(RL2),  reverse(RLl,Ll),  reverse(RL2,L2),  !. 

♦split _pair(L)  pop_global(pairiistin,[X,Y]),  cons _global(pairlistoutl,X), 

cons_global(pairlistout  2 ,  Y). 

*reverse_pairs([],[3). 

*reverse_pairs( [ [ A 3] ILL] ,[[B A]IRLL] )  reverse_pairs(LL,RLL) . 


/*  Calculates  the  fit  of  a  least-squares  2D  line  through  a  set  of  points.  */ 
*lsline(PL,[M,-l 3]) 

split_pairs(PL.XL,YL),  vprod(XLtXL,XXL), 

vprod(YL.YL.YYL),  vprod(XL,YL,XYL),surnup(XL,SXL),sumup(YL,SYL), 
sumup(XXL.SXXL),  sumup(XYL,SXYL),  length(XL,N), 

M  is  ((N*SXYL)-(SXL*SYL))/((N*SXXL)-(SXL*SXL))3  is  (SYL-(M*SXL))/N,  !. 

endptl sline( [[X 1 . Y 1  ] IBPL] ,[C  1  ,C2 ,C3] ) 
last(BPL,[X2.Y2]). 

Cl  is  (Y1  -  Y2), 

C2  is  (X2  -  XI).  C3  is  ((XI  *  Y2)  -  (X2  *  Yl)),!. 

*lslinefit(PL,[Cl,C2,C3],Fit)  > 

abs(C  1  ,AC  1  ),abs(C2,AC2),  AC1  >  AC2,!, 

reverse_pairs(PL,RPL), 

lslinefit(RPL.[C2,Cl,C3].Fit). 

*lslinefit(PL,[C  1  ,C2,C3],Fit ) 

M  is  (0-CD/C2.  B  is  (0-C3)/C2, 
split_pairs(PL.XL,YL).  vprod(XL  .XL.XXL), 

vprod(YL,YL.YYL),  vprod(XL,YL,XYL),  sumup(XL,SXL),  sumup(YL,SYL), 
sumup(XXL.SXXL),  sumup(XYL,SXYL),  length(XL.N),  sumup(YYL,SYYL), 
SqFit  is  (M*M*SXXL)  +  (2*M*B*SXL)  +  SYYL  +  (N*B*B)  +  (-2  *M*SXYL) 
+  (-2  *B*SYL), 

xsqrt(SqFit,A),Dl  is  M*M+1,  xsqrt(Dl,D), 

Fit  is  A/D,!. 


/*  Iterative  vector  processing  */ 

♦sumup(L.N)  set_global(sumuplist,L).  set_global(sumupsum,0), 
iteratelist(sumupitern.sumuplist).  sumupsum(N),  !. 

♦sumupitem  pop_global(sumuplist,X),  sumupsum(N),  retract(sumupsum(N)), 
NpX  is  N+X.  asserta(sumupsum(NpX)). 


77 


*vprod(Vl,V2,V)  set_global(vprodout,[]),  set_globa](vprodinl,Vl), 
sct_global(vprodin2,V2),  iteratelist(vproditem,vprodin  1 ), 
vprodout(RV),  rcverse(RV,V),  !  . 

♦vproditem  pop_global(vprodin  1  X) ,  pop _global(vprodin2,Y),  XY  is  X*Y, 
cons_global(vprodout,XY ). 

*vsum(Vl,V2,V)  set _global(vsumout,[]), 

set_global(vsuminl  ,V  1 ),  set_global(vsumin2,V2), 
iteratelist(vsumitem,vsuminl),  vsumout(RV),  reverse(RV,V),  !. 

♦vsumitem  pop_global(vsuminl,X),  pop _global(vsumin2,Y),  XY  is  X+Y, 
cons  .global ( vsumout  ,XY). 

*vdiff(Vl,V2,V)  set_global(vdiffout,0),  set_global(vdiffinl,Vl), 
set_global(vdiffin2,V2),  iteratelist(vdiffitem,vdiffinl ), 
vdiffout(RV),  reverse(RV,V),  !. 

♦vdiffitem  pop_global(vdiffinlX),  pop_global(vdiffin2,Y),  XY  is  X-Y, 
cons_global(vdiffout,XY). 

*sprod(Vl,K,V)  set _global(sprodout,[]),  set_globa](sprodin,Vl), 

iteratelist(sproditem(K),sprodin),  sprodout(RV),  reverse(RV,V),  !. 
♦sproditem(K)  pop_global(sprodin,X),  NX  is  X*K,  cons_globa!(sprodout,NX). 
*ssum(Vl,K,V)  set_global(ssumout,[l),  set_global(ssumin,Vl), 
iteratelist(ssumitem(K),ssumin),  ssumout(RV),  reverse(RV,V),  !. 

♦ssumitem(K)  pop _global(ssumin,X),  NX  is  X+K,  cons _global(ssumout,NX). 

/*  Management  of  global  variables  as  single-argument  facts  */ 
*set_g;lobal(Name,Value)  OLDP  =..  [Name.Oldvalue],  retract(OLDP), 

P  =..  [Name,  Value],  asserta(P),  !. 

♦set _global(Name, Value)  P  =..  [Name.Value],  asserta(P),  !. 

*cons_global(Name,I)  P=..[Name,X],  call(P),  retract(P),  NP=..[Name,[IIX]], 
asserta(NP),  !. 

*pop_global(Name, Value)  P=..[Name,[ValuelL]],  call(P),  retract(P), 
NP=..[Name,L],  asserta(NP),  !. 

*incr_global(Name)  P=..[Name,X],  call(P),  retract(P).  Xpl  is  X+l, 
NP=..[Name.Xpl],  asserta(NP),  !. 

/♦  Forward-execution  iteration,  tenninates  when  a  given  list  is  empty  ♦/ 
*iteratelist(Ipred,Lname)  repeat,  iterate 2(lpred),  P=..[Lname,[]], 
call(P),  !. 

*iterate2(Ipred)  call(Ipred),  !. 

*iterate2(Ipred). 
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module  lists. 


f*  Various  list-processing  predicates  *1 

f*  First,  the  basics  */ 

*last([X],X)  . 

*last([XIL],Y) 
last(L,Y)  . 

*append([],L,L)  . 

*append([XIL],L2,[XIL33) 
append(L,L2,L3)  . 

*reverse(L,R) 

reverse2(L,[],R)  . 

*reverse2([],L,L) 

i 

*reverse2([XIL].R,S) 
reverse2(L,[XIRJ.S)  . 

/*  Predicates  defined  from  others  */ 

*unduplicate([].[]) 

! 

*unduplicate([XIL].L2) 

member(X.L),  unduplicate(L,L2)  . 

*unduplicate(lXIL],[XIL2]) 
unduplicate(L.L2)  . 

*intersection(n.L,n)- 
*intersection([XILl J.L2,[XIL31) 

member(X,L2).  !,uitersection(Ll  .L2.L3  >. 

*intersection(lXILl J.L2.L3) 
intersection(Ll  ,L2.L3). 

reaJ_intersection(  [],L,[] ) 
real_intersection( [XIL1  ].L2,fXIL3] ) 

real_pts_member(X,L2 ),  ! ,real_intersection(L  1  ,L2,L3 ). 
real_intersection( [XIL1  ],L2,L3 ) 
real_intersection(  L 1  ,L2,L3 ). 


real_pts_member( [X 1  ,Y  1  ] ,[[X2,Y2]IL] )  XI  =:=  X2,  Y1  =:=  Y2,  !. 

reaI_pts_member(X,fYIL]) 
real_pts_member(X.L)  . 
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*member(X,[XIL])  . 

*member(X,[YILJ) 
member(X,L)  . 

*singlemember(X,[XIL]) 

!. 

*singlemember(X,[YlL]) 
singlemember(X,L) . 

module  math. 

/*  Mathematical  Formulas  not  implemented  in  M-Prolog  */ 

*xsqrt(0,0) 

*xsqrt(X,Y) 

X<1,  !,  square_bisection(Y,X,X,l ). 

*xsqrt(l,l) 

!. 

*xsqrt(X,Y) 

X>1,  RX  is  1/X,  xsqrt(RX,RY),  Y  is  1/RY. 

*square_bisection(X,Y,LO,HI) 

X  is  (LO+HI)/2,  square(X,S),  close(S,Y),  !. 

*square_bisection(X,Y,LO,HI) 

MIDPOINT  is  (LO+HI)/2,  square(MIDPOENT,S),  S<Y,  !, 
square_bisection(X,Y  >HDPOINT,HI ). 

*square_bisection(X,Y,LO,HI) 

MIDPOINT  is  (LO+HD/2.  square(MIDPOINT.S),  S>=Y.  !, 
square_bisection(X,Y,LO  .MIDPOINT ). 

*close(X.Y)  :* 

D  is  X-Y.  D  >  -1.0E-6,  D  <  1.0E-6. 

*closer(X,Y)  > 

D  is  X-Y,  D  >  -1.0E-3,  D  <  1.0E-3. 

*square(X,Y) 

Y  is  X*X. 

*expon(X,Y) 

expon2(X,l,l.l,Y). 

*expon2(X,N,S,T,S) 

T<1.0E-6,  !. 

*expon2(X,N,S.T,Y ) 
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TP1  is  X/N*T,  SP1  is  S+TP1,  NP1  is  N+l,  expon2(X,NPl,SPl,TPl,Y). 

abs(X,X)  X  >=  0,!. 

abs(X,Y)  Y  is  0  -  X,!. 

cndmod.  /*  split_merge  */ 
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